Dynamics of Limit Cycle Oscillator Subject to General Noise 
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The phase description is a powerful tool for analyzing noisy limit cycle oscillators. The method, 
however, has found only limited applications so far, because the present theory is applicable only to 
the Gaussian noise while noise in the real world often has non-Gaussian statistics. Here, we provide 
the phase reduction method for limit cycle oscillators subject to general, colored and non-Gaussian, 
noise including heavy-tailed one. We derive quantifiers like mean frequency, diffusion constant, and 
the Lyapunov exponent to confirm consistency of the results. Applying our results, we additionally 
study a resonance between the phase and noise. 
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Limit cycle oscillators effectively model various sus- 
tained oscillations in many fields of science and tech- 
nology including chemical reactions, biology, electric cir- 
cuits, and lasers [H-Q- The phase reduction method 
is a powerful analytical tool which approximates high- 
dimensional dynamics of limit cycle oscillators with sin- 
gle phase variable that characterizes timing of oscilla- 
tion [J, [B[ . Since the phase is neutrally stable, phase per- 
turbations persist in time and result in various remark- 
able phenomena where weak action leads to significant 
effects, such as those addressed in the theory of synchro- 
nization d, 0] ■ While the theory of the phase reduction 
had been developed for deterministic oscillators, recent 
studies successfully extended the theory to limit cycle os- 
cillators subject to noise am and revealed that inter- 
play between nonlinearity and noise results in fascinating 
noise-induced phenomena including frequency modula- 
tion and noise-induced synchronization 12|. |13|. 



This extended phase reduction method, however, has 
found limited applications so far, since the method is ap- 
plicable only to Gaussian noise. While the noise in the 
real world often has non- Gaussian statistics, few theo- 
ries have considered nonlinear systems subject to gen- 
eral non-Gaussian noise, which has forced people to use 
the Gaussian approximation. In particular, whether the 
phase description is still valid for oscillators subject to 
non-Gaussian noise and how quantifiers of the phase dy- 
namics should be amended remains unknown. In this 
paper, we develop the phase reduction method for limit 
cycle oscillators subject to general, colored and non- 
Gaussian noise. By correctly evaluating the influence of 
amplitude perturbations up to second order in the noise 
strength, we derive the stochastic differential equation of 
phase, which allows us to study nonlinear oscillations in 
the real world without the Gaussian approximation. To 
confirm consistency of the result, we derive closed ex- 
pressions of quantifiers of the phase dynamics such as 



mean frequency, phase diffusion constant, and the Lya- 
punov exponent. The only limitation we impose is the 
weakness of the noise. Thus, the obtained results are 
applicable even when higher order moments of the noise 
diverge as long as the second order moment is finite and 
we confirm this fact numerically. As an application of 
the results, we study a limit cycle oscillator driven by a 
phase noise with a finite correlation time and show that 
amended quantifiers precisely predict resonance between 
phase and the noise. 

We start with the case of a two-dimensional limit cycle 
oscillator and then extend our results to higher dimen- 
sions and multi-component noise. One can describe the 
evolution of the system subject to noise in terms of the 
phase jj> _§nd the amplitude deviation r from the limit 
cycle 



<t> = u + (xf((p, r)r/(t), 
f = -Xr + ag{cf),r)r){t) ; 



(1) 

(2) 



here w is the cyclic frequency of unperturbed oscillations; 
A := — (ui/2n) In A and A is the Floquet multiplier of the 
cycle, i.e., A is the average amplitude relaxation rate; 
•q(t) is a normalized noise; cr -C 1 is the noise amplitude; 
f(<p,r) and g(4>,r) are 27r-periodic in </> and represent 
sensitivity of the phase and amplitude, respectively, to 
noise. The amplitude deviation is nonuniformly scaled 
so that Eq. ([2]) is not an approximation, but uniformly 
valid over the basin of attraction of the limit cycle, as we 
rigorously show in auxiliary material [Uj . 

We use o~ as an expansion parameter; cf>(t) = 4>o(t) + 
o-<f> 1 (t)+<r 2 fo(t)+-, r(t) = an{t)+a 2 r 2 (t) + ..., f{cj>,r) = 
fo{4>) + h{4>)r + and g(<f>,r) = g (4>) + gi(^)r + ... . 
From Eqs.dTJ and ©, o (t) = wt, <j)i = /o[<£o(t)]»?(t), 
and fi = — Ari + go[4>o(t)]r)(t); the latter two formulae 
provide 
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r+oa 

n(t) = / golMt) - wr]»y(i - r) e - A ^r . (4) 
Jo 

Meanwhile, the expansion of Eq. (Q} reads 

4> = lo + vfofafo + o 2 [fU4>o)<f>iV + fMnr)] + 0(a 3 ), 

here prime denotes derivative with respect to <j>. The 
right-hand part of the latter equation except for the term 
proportional to f\ (<fi) is merely the expansion of Eq. ([T]) 
with f((f>,r) replaced by /(</>, 0). Therefore, we can keep 
the equation unexpanded with respect to <fi but add the 
correction owing to r\(t); 

j> = lo + af (cj>)n{t) + c^fMrivit) + 0(a 3 ) . 

fi(4>o) r i r ](t) is small in comparison to cr/o(0)?](t), but 
makes an average contribution of the same order (because 
(4>i ) = 0). Thus, the fluctuating part of this term is not 
principal and may be omitted; 

4> ps u + ah{<t>)r]{t) + (<r 2 fi(fo)riv(t)) ■ 
Employing expression (|4]) for r%, we obtain 

r+oo 

(fi(4>o)riV(t)) = MMt)} / 9o[Mt)-"T] G^e^dr, 

Jo 

where C(r) := (rj(t)r)(t — r)) is the noise autocorrelation 
function. Finally, the reduced phase equation up to the 
leading contributions reads 

W . (5) 

Here r is replaced with ip/ui] the corrections to <fr caused 
by replacement of <pg with <fi in the integrand are cx cr 3 
and thus negligible. Remarkably, the effect of the am- 
plitude relaxation rate A can be approximately inter- 
preted as cutting-off long-term auto-correlations of noise 
if there are some, because for large At correlation func- 
tion C(t) is suppressed by the exponential factor. Thus 
A -1 determines the maximal efficient range of noise auto- 
correlation. 

For Ornstein-Uhlenbeck noise, C(r) = 7exp(— 7|r|), 
the reduced phase equation ([5]) takes the form 



9 

(7 7 



<j> = u + afo(<t>)v(t) + —!-fi(<l>) g (4>-rp)e-^diP, 
u Jo 

which coincides with the one presented in Ref.Jjl and im- 
plies the corresponding results of Refs. [1, 0, [Tj]- While 
Ref. Q considers the case of Gaussian noise, a highly 
stable limit cycle and short noise correlation times and 
Ref. [J| is limited to the case of OU noise, the present 
theory includes their results (as special cases) and ad- 
ditionally allows dealing with non-Gaussian noise, arbi- 
trary noise auto-correlation functions (including signals 



of chaotic oscillators) and arbitrary rate of amplitude re- 
laxation. 

The procedure for deriving the reduced phase equa- 
tion suggests that this equation will provide the correct 
probability density function for (f> and mean frequency 
CI = (0) up to 0(cr 2 ); 

n = w + ^ (m) £ M<f>-4)c(£)dA 

a 2 / .,„ ,_/^__- 



+—[h(4>) J o go^-^C^je—U^J (6) 

[henceforth, (...)$ = (27r)- x /„ ...#]. The noise can ei- 
ther increase or decrease the mean frequency, depending 
on features of correlation function C(r), sensitivity func- 
tions, and the cycle stability (e.g., see Fig. [21). However, 
one should verify whether the more subtle quantities — 
the phase diffusion constant D and the leading Lyapunov 
exponent Ao — can be correctly evaluated from Eq. ([5]) . 

The principal contributions to the phase diffusion are 
readily determined from Eq. ©; indeed, 
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dr /„($ f (<f> + ujt) C(t) + 0(a 4 l7) 



4>i(t) [Eq. (J^J) ] is precisely determined by terms accounted 
in Eq. ([5]); therefore, Eq. (O is completely consistent with 
the reduced phase equation. Interestingly, up to the lead- 
ing order of accuracy the phase diffusion is not affected by 
the extra amplitude terms. Thus, for instance, the ana- 
lytical results and important conclusions of Refs. [lil. Hij 
for limit cycle oscillators subject to weak noise and de- 
layed feedback control remain correct. 

For the leading Lyapunov exponent, the situation is 
more subtle. To deal with it rigorously, we consider a 
small perturbation (a — cxq exp[/x(t)], s) to the solution 
(0(t),r(i)) of Eqs. (Ill) and We have 



a Q 



s = -As + ag' [<t>(t)]a e^r](t) + ffffi[0(t)]«j(t) 



and employ the standard multiscale method adopting 
fi(t) = m(*o,*2, ■■■), d/dt = d/dt + a 2 d/dt 2 + etc. 
After some calculations, one finds the expression for the 
leading Lyapunov exponent Ao := (fi) up to 0(a 2 ): 

^2 , /•+oo 



Ao = 
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1 UJ 
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0(a 4 ) 



= ~{f'M /oft - tf) c(£\ + 0(a 4 ), (8) 
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which is consistent with the phase equation ([5]). Note, 
in the latter equations, the amplitude degree of freedom, 
which was disregarded in previous works, impacts the in- 
stantaneous growth rate of perturbations, but averages 
out to zero. Thus, on the one hand, our results demon- 
strate the importance of amplitude degrees of freedom for 
the stability of response of a general limit cycle oscillator 
even in the limit of vanishing noise; on the other hand, 
its average impact turns out to be zero up to the leading 
order of accuracy for general noise, proving that ana- 
lytical calculations and conclusions presented in [ill EH 
are valid for real situations. Notice, the negative Lya- 
punov exponent and its decrease with increase of the 
noise strength are related to the stability of the noisy 
system response in sense that it attracts trajectories (the 
phenomenon is known as noise-induced synchronization) , 
but this does not mean that the response is regular due 
to the nonzero phase diffusion. 

All the results can be extended in a straightforward 
manner to the case of an TV-dimensional dynamical sys- 
tem subject to M-component noise; 



-p1x10' 
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0=1 



<7/3//3(</>>0) w) + 2^ 



-i x J 



< J 9i3(4>-il>,0)Cp^ ) e » dip 
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D = J2 
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fptf-MCp 



ujJ 
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df p (0~l/j,O) fif, 
7T-, Ofi 



dip 



(11) 

Here indexes noise components, j does the degrees of 
freedom transversal to the limit cycle. 

Now, we address the issue of applicability of our re- 
sults for noise with diverging higher moments. Although 
the derived expressions involve only second moments of 
the noise, one has to check that possible divergence of 
higher moments does not break the entire expansion and 
influence f2, D, and Ao in the main order. 

For this reason we performed numerical simulation of 
a Hopf oscillator subject to colored noise r](t): 

A = iA+{\/2){1-\A\ 2 )A + <jt 1 , (12) 
17 = t- 1 [-t? + *(tj) £(*)], (13) 

where A is complex, the noise acts only on Re (A) 
(Eq. ([T2")l describes, for instance, lasers with optical injec- 
tion in the limit of large density of excited states; cf. [17j). 

is Gaussian and white: (£(*)£(*')) = 25 (t - f ). We 
consider normalized noises r](t) ((f? 2 ) = 1) with three 
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FIG. 1: (Color online) Hopf oscillator (|12[) subject to different 
noises; here t v = 1 and A = 2. 

(a) : correlation function C(r) for Ornstein-Uhlenbeck noise, 
which is Gaussian, (red circles) and noises with exponential 
(blue squares) and fractional rational (green diamonds) dis- 
tributions. 

(b) -(d): the numerically calculated mean frequency (red cir- 
cles) and Lyapunov exponent (blue squares) are in good agree- 
ment with Eqs. (|14l) and (|15p (solid lines) for OU noise (b) and 
noises with exponential (c) and fractional rational (d) distri- 
butions. 



kinds of distribution V(rj): 

(1) Gaussian, V\{rf) = {2ti)- 1 / 2 exp(-r} 2 /4); 

(2) exponential, V^T]) = (1/4) exp(— \t]\/2), which has 
nonzero but still finite higher cumulants; and 

(3) fractional rational function, Vj,(rj) = 7r _1 (l + ?y 2 )~ 2 , 
for which (rj 2n ) is finite only for n = 1. 

These noises are generated with employment of si(iy) = 
1, fl 2 (Tj) = y/T/t+ \ V \/2, and s 3 ( V ) = ^/(l + V 2 )/3 in 
jEq. <(T3j). For the oscillator (p~2|) . one finds fo — — sin</>, 
fi = — fo = sin0, and go = cost/); therefore, 



n = i 

2 



D = -2A = ct 5 



sin^ (1 - e _A ^) C{ip) dip , (14) 



cosipC(ip) dip . 



(15) 



For exponential and fractional rational distributions, the 
correlation function C(t) was calculated numerically. In 
Fig. [Done can see that the analytical theory is in fairly 
good agreement with results of numerical simulation both 
for noises with all moments finite (b, c) and for one with 
infinite (rf) (d). For the latter case the analytical theory 
is practically no less accurate than for the former ones 
though, for strong noise, the mismatch between theory 
and numerics is more pronounced because of large fluc- 
tuations occurring in distributions with heavy tails. 

Another important particular opportunity yielded 
by the theory we developed is the treatment of 
the effect of the phase noise, rj(t) — s/2cos[u>ot + 
^7 / C(ti)dti}. With the noise autocorrelation function 
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FIG. 2: Hopf oscillator (|12[) subject to phase noise n(t) = 

y/2 cos[w * + /*C(*i)d*x] for cr = 0.1, 7 = 0.125, A = 0.4. 
Circles: numerical simulation, solid line: analytical theory 
[Eqs. (|14|l and (|15[)]. dashed line: analytical theory disregard- 
ing the amplitude degree of freedom. 



of freedom was disregarded (e.g., 13J), remain generally 
correct. The theory provides opportunit y fo r analytical 
investigation of the reliability of neurons 191 and consis- 
tency of lasers as well as the quality of clocks, electric 
generators, lasers, etc. for general noise and general limit 
cycle oscillators. 
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C(r) = cos(ajnr) exp(— 7|r|) one can evaluate quantifiers, 
n and Ao- In Fig.[2]the results of numerical simulation for 
the Hopf oscillator [Eq. (|12p ] subject to the phase noise 
are compared to the analytical theory. Two points are 
worth emphasizing here: (i) Now we have the phase de- 
scription for general oscillators subject to noise which is 
the representative of signals of chaotic and stochastic os- 
cillators. This is important because it provides us with a 
tool to analytically investigate the synchronizing action 
of another oscillator (either chaotic or stochastic) on the 
system under consideration in general, (ii) The ampli- 
tude degree of freedom is essential here: in the graph 
for the frequency (Fig. [2]), one can see how the analytical 
theory neglecting the amplitude perturbations (dashed 
line) is far from the real observations fairly fitted by the 
theory we have developed. The most remarkable effects 
here are observed when the characteristic noise correla- 
tion time 2-7T /ojq is commensurable with the natural oscil- 
lation period of the system, that is nonsmall, meanwhile 
the earlier studies were not able to deal with such a case. 

Summarizing, we have derived the reduced phase equa- 
tion for limit cycle oscillators subject to general non- 
Gaussian noise. The derived phase equation correctly 
provides the mean frequency, the phase diffusion constant 
and the Lyapunov exponent. Since the noise-induced 
shift of the mean frequency means the shift of the reso- 
nant frequency for entrainment by external forcing 0, S] > 
our result for mean frequency is immediately relevant for 
all investigations concerning collective phenomena in net- 
works of coupled oscillators, e.g., [H, 0, [l8j], where noise 
is unavoidably present. In particular the theory is valid 
for noise which is the representative of signals of chaotic 
and stochastic oscillators and thus may provide an ac- 
curate analytical tool to investigate their synchronizing 
action. For the Lyapunov exponent, importance of the 
amplitude degrees of freedom has been proven, though 
their average impact on the system stability vanishes in 
the leading order of accuracy. This implies that the an- 
alytical theories in earlier studies on the phase diffusion 
and the Lyapunov exponent, where the amplitude degree 
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AUXILIARY MATERIAL FOR LETTER 
"DYNAMICS OF LIMIT CYCLE OSCILLATOR 
SUBJECT TO GENERAL NOISE" 

by Denis S. Goldobin, Jun-nosuke Teramae, 
Hiroya Nakao, and G. Bard Ermentrout 

Here we provide a rigorous derivation of Eqs. (l)-(2) 
of the main article [l| (and their TV-dimensional version) 
governing evolution of a general dynamical system in the 
basin of attraction of the limit cycle for the noise-free 
case. The only restriction we impose is the differentia- 
bility of the phase flux in the basin of attraction. This 
condition allows employment of Taylor series and is ful- 
filled for a general system. 

We would like to stress that Eqs. (l)-(2) are regarded 
as a conventional paradigm for limit cycle systems, and 
Ref. Q and this auxiliary materials serve the purpose to 
confirm that this paradigm, intuitively adopted by com- 
munity, is an accurate description but not simply a model 
catching key features of the oscillatory dynamics. 

2D PHASE SPACE 

First let us consider two-dimensional case. We recall 
that the phase can be introduced as the coordinate along 
the cycle so that, it grows uniformly and increases by 2tt 
for one revolution of the system. The phase is governed 
by the equation 

4> = w, (16) 

where ui — 2ir/T is the cyclic frequency of oscillations, T 
is the period. The phase can be extended to the whole 
basin of attraction of the limit cycle so that Eq. (fT6|) holds 
valid all over the basin [Ij . Let us briefly outline the ge- 
ometric explanation for this fact. We take the state Ao 
on the phase plane (see FigOJ) an d let the system evolve 
for the time period T, the new state is Ai. In its turn, 
Ai evolves for the same time period to A2, and so forth. 
The sequence of A„ tends to the limit cycle and A^ 
belongs to it. One can connect points A and Ai by 
an arc, which can deviate from the linear segment con- 
necting these points. After each iteration for one period 
T the arc A„A„ + i turns into an arc, connecting points 
A n +i and A„+2. In such a way we end up with a curve 
running through the points Ao, Ai, A2,..., Aoo. This 
curve can be assigned the value of phase <j> at point Am 
of the limit cycle; in the literature, it is referred to as 
isochron. Obviously, such definition of phase <fi is not 
unambiguous because there are infinitely many arcs con- 
necting Ao and Ai; however it becomes unambiguous 
when one claims the curve running through A„ to be 
smooth. Possibility to construct the field of phase 4> all 
over the attraction basin is a well established fact and the 




FIG. 3: Sketch of construction of the field of (f>. 



phase field were, for instance, numerically reconstructed 
for the entire phase plane of the FitzHugh-Nagumo os- 
cillator in Q. 

Now we have to complete construction of the coordi- 
nate grid with the coordinate measuring the deviation 
from the limit cycle. It is frequently referred to as the 
amplitude. However, one has to keep in mind that it is 
rather deviation of the amplitude from the value corre- 
sponding to the limit cycle, but not the conventional am- 
plitude. Let p measures the length along the isochrones 
and p = features the position on the limit cycle. There- 
fore, 

= w, (17) 
p = F(hp) = -\(<f>)p + F 2 y>,p), (18) 

where F2 (4>, p) is a function which's Taylor series with 
respect to p starts with the term oc p 2 or higher powers. 

We want to scale the variable p: we replace variable p 
with r = h(4>, p) such that Eq. (JT5J) turns into 

r = —fir, (19) 

where fi = (2tt)^ 1 J* X((f>)d(f> is the average amplitude 
decay rate near the cycle. Now we have to reconstruct 
the function h{<j>,p) from Eqs. ([T7 | - (fT9|) . 

Let us consider isophase line <f> = 0. The phase flow 
induces mapping for p on this line; the state p(t = 0) = po 
on this line evolves after one revolution to 

00 

p(t = 2w/cj)=G(p ) =:J2GnPS; 

n=l 

the function G(po) has to be calculated from integration 
of the equation system dTTJ) and ([18]). From Eq. p^|) . 

r(t = 2tt/uj) = Ar , 

where A := exp(— 2-Kp/u). 

Matching the maps for p and r, one can find the func- 
tion ho(p) :— h{4> — 0,p). Indeed, 

r(2Tr/u) = h {p{2-K/w)) = h (G(p Q )) = h a {G(h { - 1] (r ))), 



6 



where h. 



(-i) 
o 



is the inverse function of /iq. Function 



ho(p) is a one-to-one (monotonously growing) function 
and, therefore, ^ is well-defined. On the other hand, 
r(2ir/u) = Ar . Equating two expressions for r(2n/u}) 

equality, one finds 

G(4~ 1} M) 



and applying the function hi to the both sides of the 



4- 1} (Ar ) 



(20) 



This equation can be resolved, e.g., in terms of Taylor 



series G{r) = J2^=i G n r n and 1} (r) = 



n=l 1 



n=l 



TCI — 1 



5> n A%". 



(21) 



Now one have to collect and equate terms with equal 
powers of tq. Thus, 

(1) for ro: G\a\ — a\K = 0. 

here we found an obvious claim, G\ — A, which follows 
from the fact that linearized in p evolution of small devia- 
tions from the limit cycle is determined by the multiplier 
A. Coefficient a\ remains undetermined because, in fact, 
no scale for r is imposed by our construction and we are 
free to choose oi = 1. 

(2) forr 2 : Gia 2 + G 2 a\ - a 2 A 2 = 0. 
Hence, 



0-2 



G 2 a\ (ji=i 



G, 



Gi - A 2 



(l-A)A 



(3) for rl: 



2G 2 aia 2 + G 3 af ai =i 2G 2 a 2 + G 3 



(1- A 2 )A 



(1-A 2 )A ' 



Hereby one can reconstruct the sequence of a n up to 
required order of accuracy. With the function X \p) 
evaluated one can find its inverse function ho(p). 
Let us now consider Eq. (|19p with r = h((ft,p); 

; dh dh 
P- 



K<t>, p) = <t>-^ + pf-p = h ^ p)- 

Substitution of <ft and p from Eqs. (fTTl) and (IT%1) yields 

-fih((f>,p). 



dh dh 

u 8* + F{ *' p) Tp 



One can consider this as an evolution equation 

dh p F(<f>,p) dh 

7ui = H<P,P) 7T 

o<p uj ui op 

with initial condition 

h( ( p = 0,p) = h o (p) 
calculated from Eq. (|2U| or ([2~Tj) as described. 



For small deviations p (or r), when one can neglect 
nonlinear terms in F(p) and G(p), one finds ho(p) = p 
and solution to Eq. (f22|) : 



r = h(cft, p) = pexp 



1 



(A(V) - PW 



0{P 2 )- (23) 



The particular result for small deviations, Eq. (|2"3"|) . can 
be found in Q • This result means that with an appropri- 
ate choice of the coordinates one can obtain a constant 
decay rate for amplitude deviations even when in ordi- 
nary coordinates one can observe positive instantaneous 
Lyapunov exponent (— X(cft)) meaning local divergence of 
trajectories. However, here we have shown the regular 
procedure for constructing parameterization (eft, r) such 
that, the evolution of two-dimensional dynamical system 
is accurately described by equations 



-pr, 



not only for small deviations, but all over the attraction 
basin of the limit cycle. In the presence of noise crr)(t) 
one finds 



<j> = bj + a Z^((j>, r)r)(t), 
r = -pr + a Z r (cft, r) r](t) 



(24) 
(25) 



where Z^ and Z r are sensitivity functions. This is the 
equation system (l)-(2) of the main article |l| up to no- 
tations. 



iV-D PHASE SPACE 

In higher dimensions we restrict our consideration to 
the case of small deviations because the procedure for 
consideration of nonlinearities is principally the same 
as for the two-dimensional case, but significantly more 
lengthy. The deviation from the limit cycle is now param- 
eterized by (N — l)-dimensional vector p. For linearized 
case, one can choose the point eft — and construct the 
linear mapping A; 

P{4> = 0) = A • p(<ft = 1) . 

As long as matrix A possesses only positive or complex 
eigenvalues (multipliers) Aj — exp(2irpj /lo) with eigen- 
vectors Pj, one can choose the coordinate grid (eft, f) such 
that p(t) = T.j^p^e-^, where pj(<f>) = p(t = 



(22) 4>)\p(t=o)=p r Then 

r 'j = -Pi r j +0(r 2 ). 

This equation system was assumed for derivation of the 
phase reduction equation (9) and quantifiers of the dy- 
namics (10) and (11) in the main paper. 
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In fact, our constructions correspond to the employ- 
ment of the basis of Floquet eigenvectors and develop- 
ment of this methodology for the case of nonlinear equa- 
tions. 
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